#######################################################################################
# Summary of Simulation Results: Monte Carlo Study 1 
#################################################################################

# Function to Average Across S=1000 Simulations

MeanOutFun <- function(IN_FILE, OUT_FILE){

result.array <- readRDS(paste(IN_FILE, ".rds", sep=""))
mean.array <- matrix(NA, 6, 8)
rownames(mean.array) <- rownames(result.array)
colnames(mean.array) <- c("NumRob", "TP", "FP", "TN", "FN", "CSign", "RMSE", "Corrs")
for(i in 1:6){
  for(j in 1:8){
    
    mean.array[i, j] <- mean(result.array[i, j, ], na.rm=T) 
    
  }
}

#write.table(round(mean.array, 2), file=paste("/Users/richardtraunmuller/Dropbox/PSRM_Replication/", OUT_FILE,  ".csv", sep=""), sep=",")
print(mean.array)
}

################################################################################
# Table 1: Probability of False Positive and False Negatives When Variables are Uncorrelated

p <- 1

n1 <- round(1-MeanOutFun("one_result_mat_1", "one_means_25")[c(1:2, 6, 3:5), c(2)]/p, 3)
n2 <- round(1-MeanOutFun("one_result_mat_2", "one_means_50")[c(1:2, 6, 3:5), c(2)]/p, 3)
n3 <- round(1-MeanOutFun("one_result_mat_3", "one_means_75")[c(1:2, 6, 3:5), c(2)]/p, 3)
n4 <- round(1-MeanOutFun("one_result_mat_4", "one_means_10")[c(1:2, 6, 3:5), c(2)]/p, 3)
n5 <- round(1-MeanOutFun("one_result_mat_5", "one_means_125")[c(1:2, 6, 3:5), c(2)]/p, 3)

p1 <- round(1-MeanOutFun("one_result_mat_1", "one_means_25")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p2 <- round(1-MeanOutFun("one_result_mat_2", "one_means_50")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p3 <- round(1-MeanOutFun("one_result_mat_3", "one_means_75")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p4 <- round(1-MeanOutFun("one_result_mat_4", "one_means_10")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p5 <- round(1-MeanOutFun("one_result_mat_5", "one_means_125")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)

n <- t(rbind(n1, n2, n3, n4, n5))
p <- t(rbind(p1, p2, p3, p4, p5))

tab.1.a <- rbind(n[1,], p[1,], n[2,], p[2,], n[3,], p[3,], n[4,], p[4,], n[5,], p[5,], n[6,], p[6,])

#
p <- 5

n1 <- round(1-MeanOutFun("five_result_mat_1", "five_means_25")[c(1:2, 6, 3:5), c(2)]/p, 3)
n2 <- round(1-MeanOutFun("five_result_mat_2", "five_means_50")[c(1:2, 6, 3:5), c(2)]/p, 3)
n3 <- round(1-MeanOutFun("five_result_mat_3", "five_means_75")[c(1:2, 6, 3:5), c(2)]/p, 3)
n4 <- round(1-MeanOutFun("five_result_mat_4", "five_means_10")[c(1:2, 6, 3:5), c(2)]/p, 3)
n5 <- round(1-MeanOutFun("five_result_mat_5", "five_means_125")[c(1:2, 6, 3:5), c(2)]/p, 3)

p1 <- round(1-MeanOutFun("five_result_mat_1", "five_means_25")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p2 <- round(1-MeanOutFun("five_result_mat_2", "five_means_50")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p3 <- round(1-MeanOutFun("five_result_mat_3", "five_means_75")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p4 <- round(1-MeanOutFun("five_result_mat_4", "five_means_10")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p5 <- round(1-MeanOutFun("five_result_mat_5", "five_means_125")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)

n <- t(rbind(n1, n2, n3, n4, n5))
p <- t(rbind(p1, p2, p3, p4, p5))

tab.1.b <- rbind(n[1,], p[1,], n[2,], p[2,], n[3,], p[3,], n[4,], p[4,], n[5,], p[5,], n[6,], p[6,])

#
p <- 9

n1 <- round(1-MeanOutFun("nine_result_mat_1", "nine_means_25")[c(1:2, 6, 3:5), c(2)]/p, 3)
n2 <- round(1-MeanOutFun("nine_result_mat_2", "nine_means_50")[c(1:2, 6, 3:5), c(2)]/p, 3)
n3 <- round(1-MeanOutFun("nine_result_mat_3", "nine_means_75")[c(1:2, 6, 3:5), c(2)]/p, 3)
n4 <- round(1-MeanOutFun("nine_result_mat_4", "nine_means_10")[c(1:2, 6, 3:5), c(2)]/p, 3)
n5 <- round(1-MeanOutFun("nine_result_mat_5", "nine_means_125")[c(1:2, 6, 3:5), c(2)]/p, 3)

p1 <- round(1-MeanOutFun("nine_result_mat_1", "nine_means_25")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p2 <- round(1-MeanOutFun("nine_result_mat_2", "nine_means_50")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p3 <- round(1-MeanOutFun("nine_result_mat_3", "nine_means_75")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p4 <- round(1-MeanOutFun("nine_result_mat_4", "nine_means_10")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)
p5 <- round(1-MeanOutFun("nine_result_mat_5", "nine_means_125")[c(1:2, 6, 3:5), c(4)]/(10-p), 3)

n <- t(rbind(n1, n2, n3, n4, n5))
p <- t(rbind(p1, p2, p3, p4, p5))

tab.1.c <- rbind(n[1,], p[1,], n[2,], p[2,], n[3,], p[3,], n[4,], p[4,], n[5,], p[5,], n[6,], p[6,])

table.1 <- rbind(tab.1.a, tab.1.b, tab.1.c)
rownames(table.1) <- rep(rep(c("EBA_Leamer", "EBA_sign", "EBA_90", "MA_unig", "MA_Sala", "MA_BIC"), each=2), 3)
colnames(table.1) <- c(".25", ".50", ".75", "1.00", "1.25")

write.table(table.1, file="pluemper_traunmueller_table1.csv", sep=",")

print(table.1)
##################################################################################################################################
